\chapter{马尔可夫链蒙特卡洛与Metropolis-Hastings算法}
\minitoc

\section{起点：概率建模与贝叶斯推断}

{\color{red}\textbf{[引导问题]: 我们究竟在“采样”什么？概率分布从哪里来？}}

任何蒙特卡洛算法（利用随机样本逼近积分或期望的技术）都从概率建模开始。在贝叶斯统计中（把未知量视作随机变量、用数据更新先验信念的推断框架），我们关心的是后验分布
\[
p(\theta\mid D) = \frac{p(D\mid \theta)\,p(\theta)}{p(D)},
\]
其中$\theta$是模型参数，$D$是观测数据。\textbf{似然} $p(D\mid\theta)$刻画参数对数据的解释程度，\textbf{先验} $p(\theta)$表达我们在观测之前的信念，\textbf{后验} $p(\theta\mid D)$则是融合二者后的信念更新。分母
\[
p(D) = \int p(D\mid\theta)p(\theta)\,\mathrm{d}\theta
\]
被称为\textbf{证据}或\textbf{归一化常数}。在高维参数空间中，计算这项积分往往不现实，这正是MCMC要解决的问题：\textbf{在不显式计算$\int$的情况下，从$p(\theta\mid D)$获取样本。}

{\color{red}\textbf{[思考任务]:} 当先验与似然都可评估但积分难算时，为什么我们仍然可以比较$\tilde{\pi}(\theta)=p(D\mid\theta)p(\theta)$的大小？}

{\color{red}\textbf{[参考答案]:} 我们真正希望抽样的是后验$\pi(\theta\mid D)$，它把先验知识与数据证据结合到同一个概率模型中。由于归一化常数$p(D)$对所有候选$\theta$都相同，比较$\tilde{\pi}(\theta)$的大小等价于比较后验的相对大小，因此即便无法算出积分，也能据此设计抽样与接受概率。}

\section{蒙特卡洛积分的基本原理}

{\color{red}\textbf{[引导问题]: 大数定律告诉我们“样本均值逼近期望”，那在高维空间是否依旧可靠？}}

经典蒙特卡洛积分把期望写成
\[
\mathbb{E}_\pi[f(X)] = \int f(x)\pi(x)\,\mathrm{d}x \approx \frac{1}{N}\sum_{i=1}^N f(x^{(i)}),
\]
其中$x^{(i)}\sim\pi(x)$相互独立同分布。大数定律保证当$N\to\infty$时，样本均值几乎必然收敛到真实期望，中心极限定理进一步说明误差尺度约为$\mathcal{O}(1/\sqrt{N})$。

然而，在高维问题中我们既无法均匀覆盖空间，也难以直接从$\pi$抽样。重要抽样（先从易采样的$q$生成点，再用权重校正以还原目标分布）、拒绝抽样（利用包络分布对候选点执行接受/拒绝决策）等方法依赖一个易采样的提议分布$q(x)$。一旦$q$与$\pi$差异过大，权重$w(x)=\pi(x)/q(x)$就会剧烈波动，导致估计方差爆炸。\textbf{维度灾难}意味着我们必须寻找新的方式，让样本自适应聚焦在高概率区域，而不是“盲目撒网”。

{\color{red}\textbf{[动手问题]:} 回忆一下重要抽样估计$\mathbb{E}_\pi[f(X)]$的公式，思考当$q$的支撑小于$\pi$时会发生什么。}

{\color{red}\textbf{[参考答案]:} 大数定律同样适用于高维，只要样本独立同分布，不过在实际中我们难以直接获得来自$\pi$的独立样本，因此收敛速度受限于能否有效探索高概率区域。若重要抽样的提议$q$在某些区域的密度为零，而目标$\pi$仍有正概率，则对应样本的权重$\pi/q$会发散，导致估计偏差和方差都无法控制——这正体现支撑不完全时重要抽样失效。}

\section{马尔可夫链基础：从随机游走到平稳分布}

{\color{red}\textbf{[引导问题]: 如何设计一个“无记忆”的随机过程，让它最终忘记出发点，只记得目标分布？}}

设$\{X_t\}_{t\ge 0}$是状态空间$\mathcal{X}$上的随机过程，若满足
\[
\mathbb{P}(X_{t+1}=j \mid X_t=i, X_{t-1},\ldots,X_0) = \mathbb{P}(X_{t+1}=j \mid X_t=i) = P_{ij},
\]
就称为\textbf{马尔可夫链}。矩阵$P=(P_{ij})$是转移概率矩阵，满足每一行求和为1。通过矩阵乘法可得$k$步转移概率$P^{(k)}=P^k$，这正是Chapman--Kolmogorov方程的矩阵形式。

若对任意$i,j\in\mathcal{X}$存在$k$使$P^{(k)}_{ij}>0$，链是\textbf{不可约}的；如果存在某状态$i$使得$\gcd\{k\mid P^{(k)}_{ii}>0\}=1$，则称该状态\textbf{非周期}。当链从任一状态出发的期望回返时间有限时，我们称其\textbf{正返}。对有限状态的不可约非周期链，\textbf{遍历定理}断言：无论初始分布如何，$X_t$都会收敛到唯一的\textbf{平稳分布}$\pi$，并且时间平均收敛到空间平均：
\[
\frac{1}{n}\sum_{t=1}^n f(X_t) \xrightarrow{a.s.} \sum_x f(x)\pi(x).
\]

{\color{red}\textbf{[思考问题]: 为什么“不可约+非周期”足以保证极限分布唯一？你能用地铁网络的连通性为它构造几何直觉吗？}}

在连续状态空间上，我们使用转移密度$P(y\mid x)$描述一步跳跃。只要马尔可夫链是遍历的（不可约、非周期、正返——意味着可以互通并在有限期望时间内回返），从任意初值出发的长链都可以看作来自$\pi$的相关样本。

构造MCMC时，我们常使用\textbf{细致平衡}条件：若$\pi(x)P(x,y)=\pi(y)P(y,x)$对所有$x,y$成立，则$\pi$必为平稳分布。这提供了设计转移核的直接途径，也是MH算法和HMC算法（Hamiltonian Monte Carlo，借助梯度信息模拟“动力学”跃迁）的理论基础。
\section{从疑问出发：为什么要用MCMC？}

{\color{red}\textbf{[引导问题]: 假如你手里只有一台计算器，却要评估一个100维积分，你会怎么办？}}

统计推断、贝叶斯学习乃至统计物理都有一个共同难题：我们需要从一个复杂的概率分布中抽样，而该分布只知道一个未归一化的密度函数$\tilde{\pi}(\theta)$。换句话说，我们知道$p(\theta\mid D) \propto \tilde{\pi}(\theta)$，但归一化常数$Z = \int \tilde{\pi}(\theta)\,\mathrm{d}\theta$难以计算。直接计算这个积分在高维空间中几乎不可能，传统的网格积分、梯度优化甚至重要抽样都可能因为“维度灾难”失效。

{\color{red}\textbf{[思考提示]: 如果直接抽样做不到，能不能让一个随机过程帮我们“走遍”整个空间，从而间接还原分布呢？}}

这正是马尔可夫链蒙特卡洛（Markov Chain Monte Carlo, MCMC）的基本想法：与其一次性找到独立样本，不如设计一个随机游走过程，让它在长时间的运行后，访问状态的频率与目标分布成正比。你可以想象自己在一片崎岖的山地中行走，如果你的行走规则恰到好处，那么在漫长的旅程结束时，你停留在每个高度的时间比例就反映了该高度对应的概率密度。

{\color{red}\textbf{[参考答案]:} 100维积分难以直接评估时，可以把问题交给随机过程：通过构造遍历的马尔可夫链，让样本在高概率区域多停留、在低概率区域少逗留，就能利用时间平均逼近期望。换言之，随机游走提供了间接复原目标分布的途径。}

\section{马尔可夫链视角：把采样问题转化为动力学设计}

{\color{red}\textbf{[引导问题]: 既然我们要靠随机游走帮忙，怎样的游走才能保证“不偏不倚”呢？}}

核心思想是构造一个马尔可夫链$\{X_t\}_{t\ge 0}$，其状态空间与我们想要采样的变量相同。当链的转移概率满足一定条件（不可约、非周期、正返），它会收敛到一个平稳分布$\pi(x)$。如果我们把这个平稳分布设计成目标分布，那么长时间运行后，$X_t$的经验分布就逼近$\pi(x)$。

{\color{red}\textbf{[关键概念回顾]:}
\begin{itemize}
\item \textbf{不可约性：} 从任一状态出发，都存在有限步可以到达任何其他状态。否则链会被困在子空间，永远看不到目标分布的某些部分。
\item \textbf{非周期性：} 避免链以固定节奏来回跳动而不收敛到平稳分布。
\item \textbf{细致平衡（Detailed Balance）：} $\pi(x)P(x,y)=\pi(y)P(y,x)$，它提供了一条构造平稳分布的“捷径”。
\end{itemize}}

{\color{red}\textbf{[思考问题]: 为什么细致平衡足以保证$\pi(x)$是平稳分布？}} 提示：把所有状态的概率流入和流出写出来，观察在细致平衡条件下它们如何逐项抵消。

{\color{red}\textbf{[参考答案]:} 设计马尔可夫链时关键是让转移规则既能遍历状态空间又不偏离目标分布。细致平衡逐对地平衡了$\pi(x)P(x,y)$与$\pi(y)P(y,x)$的概率流量，因此对任意状态，净流入等于净流出，整条链自然保持$\pi$不变，这就保证了$\pi$是平稳分布。}

\section{提议核与转移机制的基础知识}

{\color{red}\textbf{[引导问题]: 字母$T(y\mid x)$究竟代表什么？为什么不同的提议核会带来截然不同的链行为？}}

在MCMC中，我们常把\textbf{提议核}$T(y\mid x)$写成一个条件分布：给定当前状态$x$后，随机机制根据$T(y\mid x)$生成候选$y$。在离散情形，$T(y\mid x)$是条件概率；连续情形则是条件密度，并满足对每个$x$都有$\int T(y\mid x)\,\mathrm{d}y=1$。提议核与接受率共同决定一步转移的概率$P(x,y)=T(y\mid x)\alpha(x,y)$。

{\color{red}\textbf{[关键概念]:}}
\begin{itemize}
\item \textbf{支撑（support）匹配：} 若存在$y$使$\pi(y)>0$但$T(y\mid x)=0$，链永远无法访问该区域，遍历性被破坏。因此常要求$T$的支撑覆盖目标分布的高概率区域。
\item \textbf{对称提议（symmetric proposal）：} 当$T(y\mid x)=T(x\mid y)$时称为对称提议，例如随机游走高斯$\mathcal{N}(x,\sigma^2 I)$。此时MH接受率化简为$\min\{1,\tilde{\pi}(y)/\tilde{\pi}(x)\}$。
\item \textbf{独立提议（independence sampler）：} 若$T(y\mid x)=q(y)$与当前状态无关，就得到独立采样器。它易于分析，但若$q$与$\pi$差异大，接受率会显著下降。
\item \textbf{随状态调节的核：} 在高维问题中，常让$T$依赖局部几何，如使用协方差矩阵$\Sigma(x)$或梯度信息（Langevin/Hamiltonian提议），以提升跨越窄峡谷的能力。
\end{itemize}

{\color{red}\textbf{[实践提醒]:}}
\begin{itemize}
\item \textbf{步长尺度：} 对称随机游走中协方差$\sigma^2 I$的大小决定探索半径；过小易陷入强相关，过大则导致频繁拒绝。
\item \textbf{可调参数：} 提议核往往带有易调参数（步长、协方差、自由度等）。真实应用中，可以用试运行或自适应策略（如更新$\Sigma$或局部梯度步长）寻找合适设定。
\item \textbf{组合核：} 可以对多个子核取混合或分块更新，如$T = w T_1 + (1-w)T_2$或按参数分组逐块抽样，以兼顾局部微调与全局跳跃。
\end{itemize}

当我们理解了提议核的角色，就能更有针对性地选择、调节或自适应设计$T(y\mid x)$，让MH接受率公式中的每一项都有实际意义：$T$负责生成候选，$\tilde{\pi}$衡量目标的偏好，两者共同塑造链的混合速度。
\section{从细致平衡推导Metropolis-Hastings接受率}

{\color{red}\textbf{[引导问题]: 如果我们已经有一个提议分布$T(y\mid x)$，怎样设计“接受-拒绝”步骤才能让目标分布稳定不变？}}

设我们从当前状态$x$使用提议机制$T(y\mid x)$抽取候选$y$，接着以概率$\alpha(x,y)$接受它，否则停留在$x$。链的转移概率就是
\[
P(x,y) = T(y\mid x)\alpha(x,y),\quad x\neq y;\qquad P(x,x) = 1 - \sum_{z\neq x}T(z\mid x)\alpha(x,z).
\]

我们希望$\pi(x)$满足细致平衡：
\[
\pi(x)T(y\mid x)\alpha(x,y) = \pi(y)T(x\mid y)\alpha(y,x).
\]
由于$\alpha$的取值在$[0,1]$，一个对称的选择是
\[
\alpha(x,y) = \min\left\{1, \frac{\pi(y)T(x\mid y)}{\pi(x)T(y\mid x)}\right\},
\]
这正是Metropolis-Hastings（MH）算法的核心公式。


{\color{red}\textbf{[引导任务]: 思考为什么接受率$\alpha$只需取两边之比的最小值即可满足细致平衡。}}

{\color{red}\textbf{[参考答案]:} 通过构造满足细致平衡的转移核，可以保证链在目标分布下的概率流入与流出完全抵消，因此无记忆随机过程最终只“记得”目标分布而忘记起点。不可约性确保状态空间互联，非周期性避免链被节奏锁死，两者合起来由遍历定理保证极限分布唯一。接受率选取$\alpha(x,y)=\min\{1,\pi(y)T(x\mid y)/[\pi(x)T(y\mid x)]\}$时，若流向$y$的概率大过回流$y\to x$的那一侧就不需要再放大，从而既满足细致平衡又保持取值在$[0,1]$。}

{\color{red}\textbf{[思考问题]: 上式中分子分母的比值各自代表什么物理意义？}}

{\color{red}\textbf{[参考答案]:} 分子$\pi(y)T(x\mid y)$代表在平衡态下“从$y$回到$x$”的概率流量：$\pi(y)$刻画目标分布对状态$y$的偏好，$T(x\mid y)$补偿提议核从$y$返向$ x$的倾向。分母$\pi(x)T(y\mid x)$则对应“从$x$前往$y$”的概率流量。比较二者就相当于比较往返两侧的净流：若分子不小于分母，说明向$y$的净流量已足够大，可以把接受率截断为1；否则按比例缩放以维持细致平衡，从而保证整体稳态分布仍为$\pi$。}


\subsection*{能量视角：把概率写成玻尔兹曼因子}

统计物理提供了一个理解概率分布的经典语言：若我们把目标分布写成
\[
\pi(x) = \frac{1}{Z}\exp\bigl(-E(x)/T\bigr),
\]
其中$E(x)$称为\textbf{能量函数}（也叫势能或负对数密度），$T$是“温度”常数，$Z$是归一化常数（配分函数），则高概率区域对应低能量谷地。概率与能量的互换关系常通过定义$E(x)=-T\log \tilde{\pi}(x)$获得，温度$T$可取1（统计推断中常用无量纲温度），或在模拟退火、并行温度等算法中刻意设定不同的温度以平滑能量景观。

这个视角带来几条直觉：
\begin{itemize}
\item \textbf{能量差驱动接受度}：比较两个状态$x,y$时，$E(y)-E(x)= -T\log[\pi(y)/\pi(x)]$，能量下降意味着概率提高，反之亦然。
\item \textbf{玻尔兹曼因子}：在热力学中，状态以$\exp(-\Delta E/T)$的比例在平衡态出现；MH接受率沿用同样的逻辑，保证链在平衡时遵循目标分布。
\item \textbf{温度调节探索}：增大$T$会放大能量起伏的可接受性，使链更易跨越高能障碍；$T\to 0$则趋向贪心搜索，只保留最低能量点。
\end{itemize}

基于同样的物理类比，\textbf{模拟退火法}（Simulated Annealing）把优化问题视作在能量地形中寻找最低谷。它先在较高温度下运行接受率类似$\min\{1,\exp(-\Delta E/T)\}$的随机搜索，以容忍“上山”跳跃避开局部极小；随后逐步降低温度，减少对高能状态的接受，最终在低温极限近似收敛到全局最优。能量视角因此成为连接MCMC采样与随机优化的桥梁。

因此，把概率写成能量形式不仅连接了统计物理的直觉，也帮助我们理解后续推导中的指数项、比值与最小值运算。

{\color{red}\textbf{[进阶任务]:}} 试着把$\pi(\cdot)$写成$e^{-E(\cdot)/T}$的形式，看看接受率如何转化为能量差$\Delta E$。

{\color{red}\textbf{[参考答案]:} 将目标密度写成玻尔兹曼形式$\pi(x)=Z^{-1}\exp\bigl(-E(x)/T\bigr)$后，MH的比值变为
\[
\frac{\pi(y)T(x\mid y)}{\pi(x)T(y\mid x)} = \exp\Bigl(-\frac{E(y)-E(x)}{T}\Bigr)\times\frac{T(x\mid y)}{T(y\mid x)}.
\]
令$\Delta E = E(y)-E(x)$，拒绝/接受规则就写成
\[
\alpha(x,y) = \min\Bigl\{1, \exp\bigl(-\Delta E/T\bigr)\times\frac{T(x\mid y)}{T(y\mid x)}\Bigr\}.
\]
因此能量差$\Delta E$描述从$x$跳到$y$所需“爬坡”或“下坡”的高度：当$\Delta E<0$（能量下降）时，指数项超过1而被截断，意味着提议必然接受；当$\Delta E>0$（能量上升）时，接受概率以$\exp(-\Delta E/T)$衰减，好似玻尔兹曼因子为状态升温“付出”概率。此外，非对称提议核通过$T(x\mid y)/T(y\mid x)$修正这种偏好，保证细致平衡。若提议对称，则该比值为1，公式退化为熟悉的模拟退火接受准则$\min\{1,\exp(-\Delta E/T)\}$，清楚地展示了能量视角与MH之间的联系。}

\section{Metropolis-Hastings算法逐步拆解}

{\color{red}\textbf{[引导问题]: 若要从零实现MH采样器，你会安排哪些核心步骤，确保链既能探索又保持目标分布？}}

可以把算法拆成下面几个自然步骤：

\begin{enumerate}
\item \textbf{初始化：} 任取起始点$x^{(0)}$。它不需要是高概率区域，但决定了前期的“热身时间”。
\item \textbf{提出候选：} 给定当前状态$x^{(t)}$，从提议分布$T(y\mid x^{(t)})$抽取候选$y$。常见选择是对称分布，如高斯$y\sim\mathcal{N}(x^{(t)},\sigma^2 I)$。
\item \textbf{计算接受率：}
\[
r = \frac{\tilde{\pi}(y)T(x^{(t)}\mid y)}{\tilde{\pi}(x^{(t)})T(y\mid x^{(t)})},\qquad \alpha = \min\{1,r\}.
\]
注意这里只需要未归一化密度$\tilde{\pi}$，归一化常数自动抵消。
\item \textbf{接受或拒绝：} 生成$u\sim\mathrm{Uniform}(0,1)$。若$u\le \alpha$，令$x^{(t+1)}=y$；否则保持$x^{(t+1)}=x^{(t)}$。
\item \textbf{重复：} 重复步骤2--4，直至获得足够长的链。
\end{enumerate}

{\color{red}\textbf{[细节提示]:}
\begin{itemize}
\item 当提议分布对称时，$T(x\mid y)=T(y\mid x)$，接受率化简为$\min\{1,\tilde{\pi}(y)/\tilde{\pi}(x)\}$，这就是经典的Metropolis算法。
\item 提议分布越“跳跃”，探索越快，但拒绝率也可能上升；反之跳太小会导致缓慢“爬行”。调节步长是实践中的重要技巧。
\end{itemize}}

{\color{red}\textbf{[参考答案]:} 为了既探索又保持目标分布，最重要的是维持细致平衡：提议核$T(y\mid x)$负责产生候选，使链有机会走向新的区域；接受率$\alpha(x,y)$则用目标密度比值（必要时还要乘以$T(x\mid y)/T(y\mid x)$）来校正提议偏向，确保往返流量平衡。实际实现时，按照“初始化→提议→计算接受率→接受或拒绝→迭代”的流程执行，并根据混合情况调节步长或刷新提议协方差，就能兼顾探索效率与稳态正确性。}

\section{一个手把手的示例：目标分布为双峰的混合高斯}

{\color{red}\textbf{[引导问题]: 当目标分布有多个峰时，MH算法能否顺利在峰之间穿梭？}}

考虑一维目标分布
\[
\tilde{\pi}(x) = 0.3\,\exp\!\left(-\frac{(x+3)^2}{2}\right) + 0.7\,\exp\!\left(-\frac{(x-2)^2}{0.5}\right),
\]
它表示两个不同宽度、位置的高斯峰。我们使用对称提议$T(y\mid x)=\mathcal{N}(x,\sigma^2)$。

{\color{red}\textbf{[动手问题]:} 试着在纸上推导接受率的显式表达式，并思考$\sigma=0.1$和$\sigma=2$分别会带来怎样的行为。}

\begin{itemize}
\item 如果$\sigma=0.1$，新提议点多集中在当前点附近，链容易在某个峰附近徘徊，跨峰需要大量时间。
\item 如果$\sigma=2$，链可以迅速跨越峰值之间的“谷地”，但接受率下降；你可能看到大量连续拒绝导致链停滞。
\end{itemize}

具体来说，利用对称提议的性质可将MH接受率写成
\[
\alpha(x,y) = \min\left\{1,\frac{\tilde{\pi}(y)}{\tilde{\pi}(x)}\right\} = \min\left\{1,\frac{0.3\,\exp\!\bigl(-\tfrac{(y+3)^2}{2}\bigr)+0.7\,\exp\!\bigl(-\tfrac{(y-2)^2}{0.5}\bigr)}{0.3\,\exp\!\bigl(-\tfrac{(x+3)^2}{2}\bigr)+0.7\,\exp\!\bigl(-\tfrac{(x-2)^2}{0.5}\bigr)}\right\}.
\]
若链当前位于左峰$x\approx -3$，将提议写作$y=x+\eta$并在$\eta\sim\mathcal{N}(0,\sigma^2)$下展开，可得
\[
\log \frac{\tilde{\pi}(x+\eta)}{\tilde{\pi}(x)} \approx \eta\,\nabla\log\tilde{\pi}(x) + \tfrac{1}{2}\eta^2\,\nabla^2\log\tilde{\pi}(x).
\]
左峰附近梯度为零、海森负定，因而当$\sigma$很小时，$\alpha(x,y)\approx 1-\mathcal{O}(\sigma^2)$，说明局部微调几乎总被接受；但跨峰所需的位移大约$\delta=|2-(-3)|=5$，小方差高斯提议产生如此大跳跃的概率为$\Pr(|\eta|\ge\delta)\approx\exp(-\delta^2/(2\sigma^2))$，对$\sigma=0.1$几乎为零，从而解释了长时间滞留。

反之，当$\sigma$增大时，提议有较大概率落在另一峰邻域。例如从$x=-3$跳至$y=2$，上式中的分子、分母规模相近，因此一旦提议到位便趋于被接受；但较大的$\sigma$也会频繁将$y$抛向“谷地”区域（如$y\approx-0.5$），此处混合分布几乎完全由左侧高斯主导，可以直接计算
\[
\frac{\tilde{\pi}(y)}{\tilde{\pi}(x)} = \frac{0.3\exp\!\bigl(-\tfrac{(y+3)^2}{2}\bigr)+0.7\exp\!\bigl(-\tfrac{(y-2)^2}{0.5}\bigr)}{0.3\exp\!\bigl(-\tfrac{(x+3)^2}{2}\bigr)+0.7\exp\!\bigl(-\tfrac{(x-2)^2}{0.5}\bigr)} \approx \exp\Bigl(-\frac{(y+3)^2-(x+3)^2}{2}\Bigr) \ll 1,
\]
其中近似来自于$x$在左峰附近时第二项可忽略。带入$x=-3,y=-0.5$可得$\tilde{\pi}(y)/\tilde{\pi}(x)\approx 4.4\times10^{-2}$，说明这种跳跃几乎总被拒绝。总体平均接受率
\[
\mathbb{E}[\alpha] = \iint \min\left\{1,\frac{\tilde{\pi}(y)}{\tilde{\pi}(x)}\right\}\tilde{\pi}(x)\,\mathcal{N}(y\mid x,\sigma^2)\,\mathrm{d}y\,\mathrm{d}x
\]
随$\sigma$增大而下降，体现“走得远但被拒绝多”的折中。而跨峰概率可近似为$P(|\eta|>\delta) = 2\Phi(-\delta/\sigma)$（$\Phi$为标准正态分布函数上尾）；当$\sigma=2$时该值约$2.3\%$，足以偶尔跨峰；当$\sigma=0.5$时只有$1.5\times 10^{-5}$，几乎无法实现跨峰。

{\color{red}\textbf{[观察任务]:} 想象绘制运行轨迹：当链从左峰跳到右峰时，接受率突然变高还是变低？哪种迹象说明链尚未充分混合？}

{\color{red}\textbf{[参考答案]:} 双峰示例说明步长控制着跨峰能力：$\sigma=0.1$导致链久困一侧，$\sigma=2$则虽能跨越但容易被拒绝。实际绘制轨迹时，跨峰瞬间往往先出现一个接受率“谷”。原因在于链通常需要先穿越谷地：如上所示，当$x=-3,y=-0.5$时，$\tilde{\pi}(y)/\tilde{\pi}(x)\approx4.4\times10^{-2}$，所以接受率约$0.044$——图上表现为一次明显拒绝。仅当后续提议落在另一峰邻域（$y\approx2$）时，这个比值才恢复到$\mathcal{O}(1)$，接受率重新升高。若轨迹长期只在单个峰附近徘徊、或接受率持续偏低/偏高，说明链尚未充分混合。}

\section{实现细节、调参与常见陷阱}

{\color{red}\textbf{[引导问题]: 为什么现实中的MCMC代码常常包含``burn-in''、``thin''等看似额外的步骤？}}

\begin{itemize}
\item \textbf{Burn-in（热身期）：} 链刚启动时可能远离高概率区域，需要一段时间收敛到平稳分布。通常会丢弃前$B$个样本。
\item \textbf{Thinning（抽稀）：} 为了降低样本之间的相关性，可以每隔$k$步保留一个样本。但抽稀也会减少有效信息，是否需要取决于后续分析。
\item \textbf{有效样本量（ESS）：} 相关样本的方差比独立样本大。通过估计自相关函数可以计算ESS，帮助评估采样质量。
\item \textbf{混合（Mixing）：} 形容链在目标分布各区域之间移动的效率。慢混合意味着需要更长链条才能获得可靠估计。
\item \textbf{多链诊断：} 启动多条链，观察它们是否收敛到相似分布。Gelman-Rubin统计量$\hat{R}$是常用指标。
\end{itemize}

{\color{red}\textbf{[思考问题]: 如果两条链从不同起点出发，最终轨迹仍未重合，你觉得可能是什么原因？提议分布太窄？还是目标分布存在无法跨越的势垒？}}

{\color{red}\textbf{[参考答案]:} Burn-in帮助链摆脱劣质初值，抽稀则在相关性与样本效率之间取舍；有效样本量与混合指标共同评估质量。若多条链迟迟无法重合，通常意味着提议步长太小导致局部徘徊，或目标分布含有高势垒/远隔峰值，需放大步长、使用自适应方案或引入并行温度来改善探索。}

\section{进一步提升：哈密顿蒙特卡洛（HMC）}

{\color{red}\textbf{[引导问题]: 如果随机游走老是“走两步退一步”，能否借助物理规律，让采样者像滑翔一样快速穿越高维空间？}}

MH算法的效率常常受制于随机游走的低效探索：步长太大被拒绝，步长太小又爬行缓慢。哈密顿蒙特卡洛（Hamiltonian Monte Carlo, HMC）借用了经典力学里的哈密顿动力学，引入“动量”变量，让采样轨迹沿着势能形状平滑滑行，从而在高维空间中实现远距离、低拒绝率的跳跃。

{\color{red}\textbf{[关键构想]:}
\begin{itemize}
\item 为参数$\theta$配备同维度的辅助动量$p$，定义哈密顿量$H(\theta,p) = U(\theta) + K(p)$，其中$U(\theta) = -\log \tilde{\pi}(\theta)$是势能，$K(p) = \tfrac{1}{2} p^\top M^{-1} p$是动能，$M$是质量矩阵。
\item 组合分布$\pi(\theta,p) \propto e^{-H(\theta,p)}$的边缘仍为目标分布$\pi(\theta)$。
\item 通过哈密顿方程
\[\frac{\mathrm{d}\theta}{\mathrm{d}t} = M^{-1}p,\qquad \frac{\mathrm{d}p}{\mathrm{d}t} = -\nabla U(\theta)\]
模拟一段时间，就能在相空间中沿着能量等势面快速移动。
\end{itemize}}

{\color{red}\textbf{[思考问题]: 为什么哈密顿动力学能保持$H(\theta,p)$不变？}} 提示：回顾保守力场的能量守恒，以及辛结构（symplectic structure）带来的体积保持。

\subsection*{Leapfrog数值积分：保持能量的秘密武器}

实际计算中无法解析求解哈密顿方程，需要数值积分。HMC使用“Leapfrog”（蛙跳）积分器，它具备时间可逆和辛结构，因而能较好地保持能量：

\begin{enumerate}
\item 半步更新动量：$p_{t+\frac{1}{2}} = p_t - \tfrac{\varepsilon}{2} \nabla_{\theta} U(\theta_t)$。
\item 一步更新位置：$\theta_{t+1} = \theta_t + \varepsilon M^{-1} p_{t+\frac{1}{2}}$。
\item 半步更新动量：$p_{t+1} = p_{t+\frac{1}{2}} - \tfrac{\varepsilon}{2} \nabla_{\theta} U(\theta_{t+1})$。
\end{enumerate}

重复上述步骤$L$次，就得到一条长度$L\varepsilon$的拟动力学轨迹，候选状态为$(\theta',p')$。

{\color{red}\textbf{[引导问题]: Leapfrog为什么要“半步-整步-半步”交替？}} 提示：试着比较Euler积分与Leapfrog，看看哪个满足时间可逆性。

\subsection*{HMC采样流程}

完整的HMC迭代如下：

\begin{enumerate}
\item \textbf{动量重采样：} 从$p^{(t)} \sim \mathcal{N}(0, M)$独立抽样，保证动量的边缘分布正确。
\item \textbf{动力学演化：} 以$(\theta^{(t)}, p^{(t)})$为初值，使用Leapfrog积分演化$L$步，得到$(\theta', p')$。
\item \textbf{Metropolis校正：} 计算
\[
\alpha = \min\left\{1, \exp\bigl(-H(\theta',p') + H(\theta^{(t)},p^{(t)})\bigr)\right\}.
\]
由于Leapfrog近似能量守恒，上式通常接近1。
\item \textbf{接受或拒绝：} 按概率$\alpha$接受候选；若拒绝，则保持旧状态，并把动量取负（维持时间可逆性）。
\end{enumerate}

{\color{red}\textbf{[思考问题]: 如果Leapfrog步长$\varepsilon$太大、迭代步数$L$太多，会对接受率产生什么影响？}}

\subsection*{调参与实践经验}

\begin{itemize}
\item \textbf{步长$\varepsilon$：} 太大导致能量漂移严重、拒绝率升高；太小则需更多步数才能跨越长距离。常使用双平均（dual averaging，一种基于镜像下降思想的自适应步长调整法）自适应调节。
\item \textbf{轨迹长度$L\varepsilon$：} 控制每次跳跃的“飞行距离”。$L$过小退化为局部随机游走，过大则绕圈浪费计算。
\item \textbf{质量矩阵$M$：} 若参数尺度差异大，自适应$M$有助于“等速”探索。实践中常用对角或低秩近似。
\item \textbf{No-U-Turn Sampler (NUTS)：} 自动选择轨迹长度，避免手动调$L$，通过检测“折返”迹象自适应停止扩展；这是Stan等软件的默认实现。
\item \textbf{梯度计算：} HMC依赖$\nabla U(\theta)$。自动微分框架（如PyTorch、JAX，可自动构建梯度）使得实现复杂模型的HMC成为可能。
\end{itemize}

{\color{red}\textbf{[引导问题]: 面对高度相关的目标分布，为什么合适的质量矩阵能显著提升HMC效率？}}

进一步优化可引入\textbf{Riemannian预条件}：与其使用固定质量矩阵$M$，可以让它随参数位置调整，如采用Fisher信息矩阵$G(\theta)$的逆矩阵，使动量采样分布$\mathcal{N}(0, G(\theta)^{-1})$更贴合目标的局部几何。Riemannian Manifold HMC（RMHMC）沿着黎曼度量的测地线推进轨迹，理论上能在高度弯曲的后验中维持接近正交的坐标系，大幅减小相关性。实践中需解隐式方程以更新动量，计算量较大，但在强耦合参数问题上能显著提升有效样本量。

另一种常见策略是\textbf{分块更新}：将参数拆成多个子向量$\theta=(\theta^{(1)},\ldots,\theta^{(B)})$，对某些块使用HMC（或更简单的MH），对其余块采用解析抽样、Gibbs更新或独立提议。分块允许针对每个子系统定制最佳的步长、质量矩阵或预条件器，既降低了单次梯度计算的维度，又可结合领域知识（如层次模型中的局部/全局参数分离）。这种混合式MCMC常与HMC联合使用，例如对高维连续块跑短轨迹，对离散或条件分布简单的块直接采样，从而在复杂模型中兼顾效率与实现便利性。

\subsection*{与MH的对比}

{\color{red}\textbf{[思考任务]:} 用下表总结两种算法的区别，并思考何时选择HMC。}

\begin{center}
\begin{tabular}{p{0.28\textwidth}p{0.32\textwidth}p{0.32\textwidth}}
\toprule
 & \textbf{MH随机游走} & \textbf{HMC} \\
\midrule
提议机制 & 依赖本地随机扰动 & 沿梯度驱动的动力学轨迹 \\
所需信息 & 目标密度比值 & 目标密度\emph{及其梯度} \\
典型步长 & 手动调节、常较短 & 可远距离跳跃 \\
接受率 & 易受步长影响 & 通常接近1 \\
适用场景 & 梯度难算、低维 & 梯度可得、高维连续空间 \\
实现复杂度 & 低 & 中等，需要数值积分 \\
\bottomrule
\end{tabular}
\end{center}

{\color{red}\textbf{[开放问题]: 如果目标分布存在多个远隔的能量井，HMC是否总能跨越？你会考虑哪些改进（如温度扩散、多模态HMC）？}}

更严格的比较可以从随机游走与哈密顿动力学的扩散性质入手：对维度为$d$的各向同性高斯$\pi(\theta)\propto\exp(-\|\theta\|^2/2)$，Metropolis随机游走以方差$\sigma^2$步长更新，其位置增量满足$\mathbb{E}[\|\theta_{t+1}-\theta_t\|^2]\approx\sigma^2$，而为了保持固定接受率，最佳步长满足$\sigma^2=\mathcal{O}(1/d)$（Roberts等，1997），导致每迭代探索半径按$1/\sqrt{d}$衰减。相比之下，HMC在一次轨迹上累计位移近似为$\varepsilon L \|M^{-1}p\|$：若对动量$p\sim\mathcal{N}(0,M)$，则$\mathbb{E}[\|M^{-1}p\|^2]=\mathrm{tr}(M^{-1})$，在合理选择$M$（例如与后验协方差同阶）时尺度为$\mathcal{O}(d)$。因此只要$\varepsilon$和$L$不随$d$急剧缩放，HMC的探索距离可保持在$\mathcal{O}(1)$量级，远优于随机游走的$\mathcal{O}(1/\sqrt{d})$。

此外，接受率也体现了轨迹保守的优势：Leapfrog积分的能量误差$\Delta H = H(\theta',p')-H(\theta,p)$通常为$\mathcal{O}(\varepsilon^3L)$，因此Metropolis校正的拒绝概率$1-\alpha$约为$\mathbb{E}[\max(0,1-e^{-\Delta H})] \approx \mathcal{O}(\varepsilon^3L)$，远小于随机游走中步长过大导致的指数级拒绝。综合两者，HMC在高维连续空间能以更高的接受率实现大尺度跃迁，收敛时间常随维度线性增长，而随机游走则呈二次甚至更差的增长。

{\color{red}\textbf{[参考答案]:} HMC通过引入动量让链沿哈密顿动力学“滑翔”，因此能同时兼顾远距离移动与较高接受率；能量守恒源于哈密顿方程的辛结构，Leapfrog的半步-整步-半步更新正是为了保持时间可逆与体积守恒。步长过大或步数过多都会破坏能量守恒、降低接受率；合理的质量矩阵能沿相关方向伸缩步长，提高探索效率。与MH对比时要关注梯度需求、典型步长、接受率等维度。若目标含多能量井，可结合并行退火、NUTS加温、重启或多模态HMC等策略，帮助链跨越势垒。上述标度分析说明，在高维Gaussian或光滑目标下，Metropolis随机游走的最优步长缩放为$\mathcal{O}(1/\sqrt{d})$而HMC可保持$\mathcal{O}(1)$的有效移动距离，这正是其在现代贝叶斯推断中广受青睐的数学根基。}

\section{条件抽样策略：Gibbs采样}

{\color{red}\textbf{[引导问题]: 当条件分布容易处理、但联合分布复杂时，能否绕开接受-拒绝步骤，直接逐块采样？}}

Gibbs采样是MH的一个特殊情形：若我们能写出目标分布$\pi(\theta)$的各维条件分布$\pi(\theta_j\mid \theta_{-j})$，便可依次从这些条件分布中抽样而无需计算接受率。具体算法如下，设参数分成$B$个块$\theta=(\theta^{(1)},\ldots,\theta^{(B)})$：
\begin{enumerate}
\item 选择更新顺序（常见为循环或随机扫描）。
\item 对于第$b$个块，从条件分布
\[
\theta^{(b)}_{t+1} \sim \pi\bigl(\theta^{(b)}\mid \theta^{(1)}_{t+1},\ldots,\theta^{(b-1)}_{t+1},\theta^{(b+1)}_{t},\ldots,\theta^{(B)}_{t}\bigr)
\]
直接采样。
\item 重复直至所有块更新完毕，从而得到下一个联合样本。
\end{enumerate}

由于条件抽样保持联合分布不变（可视为MH的提议始终被接受，接受率为1），Gibbs自然满足细致平衡，且不会额外引入被拒绝的步骤。这与前文的“分块更新”相呼应——若某些块的条件分布可解析抽样，就用Gibbs；其余块可继续使用随机游走MH或HMC，实现混合式MCMC。

从能量视角看，条件分布对应沿坐标轴方向搜索局部能量最优点；当目标分布呈现窄而对齐的峡谷（如高斯共轭模型），Gibbs能以低成本在各个条件子空间内穿梭，从而快速降低相邻样本间的相关性。在高维高斯联合分布中可定量比较：若协方差矩阵为$\Sigma$，随机游走的最优步长随谱半径缩放$\sigma^2\propto \lambda_\text{min}(\Sigma)$，而逐坐标Gibbs的自相关系数等于对应的偏相关系数$\rho_{jk\mid -jk}$，往往远小于1，体现其对线性相关的天然适配性。

实践中：
\begin{itemize}
\item 在层次贝叶斯模型中，局部参数的条件分布多为标准族分布（正态、Gamma、Dirichlet等），可直接Gibbs抽样；全局超参数若无解析条件分布，则用HMC或自适应MH更新。
\item 数据增强或潜变量模型（如隐变量高斯过程、潜在Dirichlet分布）也常借助Gibbs，通过引入辅助变量令条件分布恢复到熟悉的形式。
\item Gibbs属于更广泛的\textbf{坐标变换/预条件}思路：在适合的变量组合下，每个子步都是“低能”移动，从而与之前讨论的Riemannian预条件和分块更新自然衔接。
\end{itemize}

{\color{red}\textbf{[思考问题]: 在一个二维目标分布中，若条件分布均为正态且方差不依赖另一维，证明Gibbs采样的转移核等价于MH使用独立正态提议且接受率恒为1。}}

{\color{red}\textbf{[参考答案]:} 对二维情形$\theta=(\theta_1,\theta_2)$，Gibbs每一步按顺序抽样
\[
\theta_1' \sim \mathcal{N}(m_1(\theta_2), v_1),\qquad \theta_2' \sim \mathcal{N}(m_2(\theta_1'), v_2),
\]
构成转移密度$P((\theta_1,\theta_2),(\theta_1',\theta_2')) = \mathcal{N}(\theta_1'\mid m_1(\theta_2),v_1)\mathcal{N}(\theta_2'\mid m_2(\theta_1'),v_2)$。若用MH一步同时提议$(\theta_1',\theta_2')$并取$T((\theta_1,\theta_2),(\theta_1',\theta_2')) = P((\theta_1,\theta_2),(\theta_1',\theta_2'))$，代入接受率公式
\[
\alpha = \min\left\{1,\frac{\pi(\theta_1',\theta_2')\,T((\theta_1',\theta_2'),(\theta_1,\theta_2))}{\pi(\theta_1,\theta_2)\,T((\theta_1,\theta_2),(\theta_1',\theta_2'))}\right\}=1,
\]
因此Gibbs可视为MH在特殊提议下的退化形式。这个观察也解释了为什么Gibbs和MH能无缝混合使用：它们共享同一细致平衡架构。}

\section{群体智能视角：遗传算法与粒子群的MCMC关联}

{\color{red}\textbf{[引导问题]:} 遗传算法（GA）与粒子群优化（PSO）都是群体智能方法，主要用于最优化。能否从MCMC的角度理解它们的搜索机制，并据此评估它们与统计采样的差异？}

\subsection*{遗传算法：选择-变异的马尔可夫化}

遗传算法维护一组个体$\mathcal{P}_t = \{x^{(1)}_t,\ldots,x^{(m)}_t\}$，每个个体以适应度$F(x)$评估，经历选择、交叉、变异生成下一代。若把整个种群视作马尔可夫状态$Z_t=(x^{(1)}_t,\ldots,x^{(m)}_t)$，则GA的演化可写成转移核$P(Z_{t+1}\mid Z_t)$。在许多实现中：
\begin{itemize}
\item \textbf{选择}阶段按适应度归一化概率抽取父代，相当于对状态施加与$F(x)$成比例的再权重；
\item \textbf{交叉/变异}提供随机扰动，可看成提议核$T(Z'\mid Z)$；
\item 若引入\textbf{模拟退火式温度$T$}，并以概率$\min\{1,\exp(-(E(Z')-E(Z))/T)\}$接受新种群，则GA可严格满足细致平衡，其中$E(Z)=-\sum_i \log F(x^{(i)})$。
\end{itemize}
标准GA通常采用“贪心”选择——直接替换成最优个体，从而偏离MCMC的可逆性；但通过引入Metropolis接受步骤或Boltzmann选择$\propto \exp(F(x)/T)$，即可把GA嵌入MCMC框架。此类混合算法有时称为\textbf{Metropolis-Hastings GA}，在多峰目标下能兼顾探索与采样。

\subsection*{粒子群优化：速度更新与随机游走}

粒子群优化维护带速度的粒子$(x_t^{(i)},v_t^{(i)})$，更新规则
\[
v_{t+1}^{(i)} = \omega v_t^{(i)} + c_1 r_1 (p^{(i)} - x^{(i)}_t) + c_2 r_2 (g - x^{(i)}_t),\qquad x^{(i)}_{t+1} = x^{(i)}_t + v_{t+1}^{(i)},
\]
其中$p^{(i)}$是粒子历史最优，$g$是全局最优，$r_1,r_2\sim U(0,1)$. 若把整个粒子集合作为状态，显然转移依赖于当前群体的全局最优$g$，因此PSO在标准形式下并非Markovian（因为$p^{(i)},g$累积过去信息）。然而，在扩展状态向量中包含$p^{(i)},g$，就能恢复马尔可夫性。此时PSO相当于带“记忆”的随机游走，驱动力指向历史最优。

与MCMC比较差异：
\begin{itemize}
\item PSO缺少可逆性约束，速度更新是确定性+噪声的混合，难以对应某个平稳分布；
\item 通过引入温度衰减、随机重置、或在更新后以Metropolis准则接受新位置，可将PSO改造为探索指定能量函数的采样器（称为\textbf{Annealed PSO}）。
\end{itemize}

\subsection*{联系与互补}

从“能量视角”来看，GA与PSO都试图在能量面$E(x)=-\log F(x)$上寻找低能量区域，其搜索行为可视作带有偏置的随机游走。与MCMC相比，它们：
\begin{itemize}
\item 使用\textbf{群体}而非单链，依靠并行探索克服局部极值；
\item 通常缺乏细致平衡与显式平稳分布，因此更像“启发式优化”而非“精确采样”；
\item 可通过加入温度、Metropolis接受步骤或对适应度归一化，使其向MCMC靠拢，进而获得采样意义（例如跨越能量障碍、估计后验尾部）。
\end{itemize}
在复杂推断中，常见的混合策略是：先用PSO或GA定位高概率区域，再以MH/HMC在该区域内做精细采样；或使用GA/PSO提出多起点，为并行MCMC提供良好初值。这些组合体现了群体智能与MCMC的互补关系。

{\color{red}\textbf{[思考任务]:} 设目标能量$E(x)$已知，尝试设计一个“Metropolis-PSO”流程：粒子依据标准PSO更新后，以$\min\{1,\exp(-(E(x')-E(x))/T)\}$的概率接受新位置。分析该算法是否满足细致平衡，并讨论温度$T\to 0$与$T\to \infty$时的行为。}

{\color{red}\textbf{[参考答案]:} 设PSO更新提供一个提议核$T(x'\mid x)$，若我们在提议后以Metropolis接受概率$\alpha(x,x')=\min\{1,\exp(-(E(x')-E(x))/T)\}$筛选新位置，则整体转移概率为$P(x\to x')=T(x'\mid x)\alpha(x,x')$。为了细致平衡，需要存在某个平稳分布$\pi(x)\propto\exp(-E(x)/T)$使得
\[
\pi(x)T(x'\mid x)\alpha(x,x') = \pi(x')T(x\mid x')\alpha(x',x).
\]
若PSO提议核可视为对称（例如随机系数$r_1,r_2$服从相同分布，且惯性项满足$\omega=\omega^{-1}$等极端设定），则上式成立；但标准PSO包含指向全局最优$g$的偏置，导致$T(x'\mid x)\neq T(x\mid x')$，因此一般不能保证严格细致平衡。实践中可通过对称化提议（双向更新、随机选择$p,g$）或在扩展状态空间$(x,p,g)$上建立对称核，使Metropolis-PSO成为合法的MCMC算法。

温度极限：
\begin{itemize}
\item $T\to 0$时，接受率趋向$\mathbf{1}_{\{E(x') \le E(x)\}}$，算法退化为贪心局部搜索，仅当能量下降才移动，容易陷入局部最优。
\item $T\to \infty$时，$\alpha\approx 1$，等价于纯PSO，无能量筛选，保留高探索性但失去采样意义。
\end{itemize}
因此合适的温度调度（如退火策略）或对称化提议对于兼顾探索与收敛至目标分布至关重要。}

\section{贝叶斯案例：洛特卡-沃尔泰拉参数推断}

{\color{red}\textbf{[引导问题]: 当模型不可解析、积分无法显式计算时，我们如何把MCMC落地到真实科学问题上？}}

参考经典的猞猁-雪兔生态系统数据（1900--1920年毛皮交易记录），我们考虑耦合微分方程
\[
\frac{\mathrm{d}x}{\mathrm{d}t} = \alpha x - \beta xy,\qquad \frac{\mathrm{d}y}{\mathrm{d}t} = -\gamma y + \delta xy,
\]
用参数$\theta=(\alpha,\beta,\gamma,\delta)$描述捕食-猎物动力学。该模型没有简单解析解，需数值积分预测$(x(t),y(t))$。

\subsection*{推断流程拆解}

\begin{itemize}
\item \textbf{似然设计}：假设观测与模型预测之间的误差服从对数正态或高斯噪声，可写出$p(D\mid\theta)$。
\item \textbf{先验设定}：在缺乏强先验知识时，可以为每个参数选择宽泛的正态或均匀分布，确保支撑覆盖生物学上合理范围。
\item \textbf{目标分布}：后验$p(\theta\mid D)\propto p(D\mid\theta)p(\theta)$，只需计算到常数因子。
\item \textbf{提议机制}：采用对称的多元高斯随机游走$T(\theta'\mid\theta)=\mathcal{N}(\theta,\Sigma)$，协方差$\Sigma$控制步长。
\item \textbf{算法实现}：按照MH步骤迭代，针对每个提议点需一次数值解ODE，计算模型-数据残差。
\end{itemize}

具体实现时，可将连续时间方程离散到观测时刻$t_1,\ldots,t_m$：给定参数$\theta$与初始状态$(x_0,y_0)$，使用Runge--Kutta等稳定的显式积分器求得预测轨迹$(\hat{x}_{\theta}(t_i),\hat{y}_{\theta}(t_i))$。若假定观测噪声为对数正态，则对每个时间点贡献的对数似然为
\[
\log p(D_i\mid \theta) = -\frac{1}{2}\left(\frac{\log x_i - \log \hat{x}_{\theta}(t_i)}{\sigma_x}\right)^2 -\frac{1}{2}\left(\frac{\log y_i - \log \hat{y}_{\theta}(t_i)}{\sigma_y}\right)^2 + C,
\]
其中$x_i,y_i$为观测数据，$\sigma_x,\sigma_y$是噪声尺度。整条轨迹的对数似然即这些项的求和。若采用高斯噪声，只需将上式中的对数替换为直接差值即可。该步骤凸显了“每次提议都要完整跑一遍模型”的计算开销，也是实践中需要关注数值稳定性的原因。

为了加速收敛，可以在若干迭代后根据当前样本方差更新协方差矩阵$\Sigma$（自适应随机游走），或对参数进行重新缩放，让ODE的灵敏度在各维度更加均衡，从而避免链沿某些方向缓慢蠕动。

{\color{red}\textbf{[思考任务]:} 观察参数轨迹图时，什么迹象表明链还未“忘记”初值？残差是否“漂移”提示提议协方差需要调节？}

\subsection*{结果解读}

\begin{itemize}
\item \textbf{轨迹图}确认链围绕高概率区域震荡，避免呈现单调趋势。
\item \textbf{后验直方图}展示参数不确定度；峰值即最大后验估计，宽度反映可信区间。
\item \textbf{预测包络}：从后验采样多组参数，绘制仿真轨迹与真实数据对照，观察多数数据点是否落入置信带。
\end{itemize}

实际分析中，可结合以下指标判断采样质量：
\begin{itemize}
\item 计算残差序列$\epsilon_i = x_i-\hat{x}_{\theta^{(s)}}(t_i)$（对$y$亦同），检查其均值是否近零、方差是否稳定，若残差呈系统正负偏移，说明似然或噪声假设需要调整。
\item 在参数空间绘制成对散点图（pair plot），观察是否出现强烈线性关联或多模态结构；若存在狭长峡谷，可尝试对参数取对数或重新定义，使后验等高线趋于椭圆形。
\item 将若干后验样本带入ODE求解，计算与真实数据的均方误差(MSE)或对数预测密度，量化模型在预测层面的解释力。
\end{itemize}

通过该案例，读者可把抽象的“马尔可夫链”具体化为：一段段与数值ODE求解交织的运算流程，亲眼见证MCMC如何在复杂模型中提供参数后验估计。

{\color{red}\textbf{[参考答案]:} 该案例说明在无法解析求解的系统里，MCMC仍可通过数值解ODE评估似然，从而得到参数后验：若轨迹图仍沿单调趋势爬行，说明初值影响未消；若残差呈系统性偏移或协方差方向上产生“棘齿”，则应重新调整提议协方差或步长。模型预测包络覆盖观测点表明推断结果能够解释数据。}

\section{诊断与评估：相关性、收敛与有效样本}

{\color{red}\textbf{[引导问题]: 生成了长链之后，如何判断结果是否可信？}}

\subsection*{蒙特卡洛误差与有效样本量}

相关样本导致估计方差增大。定义自相关$\rho_k = \mathrm{Corr}(f(X_t), f(X_{t+k}))$后，可计算\textbf{有效样本量}
\[
N_{\text{eff}} = \frac{N}{1 + 2\sum_{k=1}^{\infty}\rho_k},
\]
并估计\textbf{蒙特卡洛误差}（MC error）：$\mathrm{SE}(\bar{f}) \approx \sqrt{\widehat{\mathrm{Var}}(f)/N_{\text{eff}}}$。经验规则要求MC error小于后验标准差的5\%。

在实作中，$\rho_k$通常通过样本自协方差估计：对目标序列$\{f(X_t)\}$先减去样本均值，再计算
\[
\hat{\rho}_k = \frac{\sum_{t=1}^{N-k} (f(X_t)-\bar{f})(f(X_{t+k})-\bar{f})}{\sum_{t=1}^{N} (f(X_t)-\bar{f})^2}.
\]
由于远距离滞后往往噪声大，实际求和时会在$\hat{\rho}_k$首次出现负值后截断，或采用Geyer的初始正序列法压缩波动。若$N_{\text{eff}}$远小于样本长度，说明链高度相关，需要重新调节提议或考虑抽稀。

{\color{red}\textbf{[思考任务]:} 当自相关缓慢衰减时，缩短抽样间隔能否显著提升$N_{\text{eff}}$？还是应该改进提议分布？}

\subsection*{诊断工具箱}

\begin{itemize}
\item \textbf{轨迹图}：检查是否出现趋势或长时间卡在某一区域。理想图像像厚实的“毛毛虫”，如果呈现缓慢爬坡或大尺度漂移，说明样本尚未脱离初值。
\item \textbf{自相关函数图}：评估链的混合速度，确认是否需要抽稀或调参。ACF若在十几个滞后内迅速衰减，则相关性可接受；若拖尾缓慢，应考虑加大步长或使用更智能的提议核。
\item \textbf{Gelman--Rubin BGR统计量}：并行运行多条链，计算链间方差$B$与链内方差$W$，指标
\[
\hat{R} = \sqrt{\frac{n-1}{n} + \frac{B}{nW}}
\]
接近1意味着链间差异与链内波动相当，提示收敛。这里$n$是每条链的后期样本数；若$\hat{R}>1.1$，说明不同初值的链仍在探索不同区域，需要延长运行或改进提议。
\item \textbf{分位数对比}：估计后验区间（如2.5\%和97.5\%）并随迭代监控稳定性。将迭代划分为多个窗口，对比每段分位点是否一致，可揭示潜伏的漂移或多模态转换。
\item \textbf{能量或对数后验序列}：对能量型模型，可直接绘制$-\log \tilde{\pi}(X_t)$，若该序列展示周期性或长时间停滞，也暗示混合不佳。
\end{itemize}

在诊断时，通常遵循“多指标交叉验证”的原则：单个指标可能被偶然噪声欺骗，例如轨迹图看似稳定但$\hat{R}>1.1$仍提示链间差异，因此建议同时检查$N_{\text{eff}}$、$\hat{R}$、ACF、分位数稳定性以及领域约束。

{\color{red}\textbf{[实践建议]:}}
\begin{itemize}
\item 从分散初值启动多条链，避免都落入同一局部模式。
\item 不依赖单一诊断指标，综合轨迹图、$\hat{R}$、$N_{\text{eff}}$与领域知识判断。
\item 识别多峰分布时，可结合温度扩散（parallel tempering，在多个温度下同时运行链并偶尔交换状态）或多链交换策略改进探索。
\end{itemize}

{\color{red}\textbf{[开放问题]:} 收敛诊断只能说明“没有明显问题”，如何在高维多模态情形下设计更可靠的验收标准？}

{\color{red}\textbf{[参考答案]:} 长链是否可信取决于方差评估与多种诊断的交叉印证：有效样本量把相关性折算成独立样本数，MC error据此评估估计精度；若自相关拖尾，可以增大步长、使用对角预条件或引入独立提议。更可靠的收敛判断往往需要多链统计（如$\hat{R}$）、分块分析、子序列比较与领域约束联合使用，在多模态时还可通过并行温度、桥采样（bridge sampling，利用两条分布之间的“桥”估计归一化常数并检查探索覆盖度）或能量诊断（energy diagnostic，通过能量分布和能量差检查链是否充分穿越势垒）确认所有峰值都被访问。}

\section{回顾与进一步思考}

{\color{red}\textbf{[复盘任务]:} 用自己的语言回答以下问题，看看你是否真正理解了MH算法。}

\begin{enumerate}
\item 为什么MH算法可以只用未归一化密度$\tilde{\pi}(x)$？
\item 细致平衡是否是MH的必要条件？如果我们放宽到全局平衡会发生什么？
\item 当提议分布不对称时，接受率中的$\frac{T(x\mid y)}{T(y\mid x)}$扮演什么角色？
\item 面对强相关、维度极高的参数空间，我们可以在MH的基础上做哪些改进？（提示：考虑自适应MH、Hamiltonian Monte Carlo等）。
\end{enumerate}

{\color{red}\textbf{[开放问题]:} 你能否设计出一个在高维空间表现良好的提议分布？它需要满足哪些几何性质，才能既保持高接受率，又能快速穿梭于不同区域？}

通过这些问题的反思，你会更清楚自己对MH算法理解的深度。下一步，你可以尝试实现一个简洁的MH采样器，亲手观察参数调节如何影响混合速度和收敛效果——亲身实践往往是加深直觉的最好办法。

{\color{red}\textbf{[参考答案]:} (1) MH只依赖未归一化密度是因为归一化常数在接受率中相互抵消；(2) 细致平衡是充分条件，但不是唯一要求，只要满足全局平衡且保证遍历性，链仍会收敛，只是构造更复杂；(3) 不对称提议时$T(x\mid y)/T(y\mid x)$负责校正提议偏好，使整体保持目标分布；(4) 处理高维或相关变量时，可采用自适应MH、HMC、Riemannian预条件或分块更新。高维提议分布应贴合后验等高线（例如协方差与后验海森矩阵一致），在长薄峡谷方向放大步长、在陡峭方向缩短，以兼顾高接受率与快速穿梭。}
